Graphene-based torsional resonator from molecular dynamics simulation 
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Molecular dynamics simulations are performed to study graphene-based torsional mechanical 
resonators. The quality factor is calculated by Qf = t<jr/27r, where the frequency oj and life time r 
are obtained from the correlation function of the normal mode coordinate. Our simulations reveal the 
radius-dependence of the quality factor as Qf ~ 2628/(22ii~^ -I- 0.004i?^), which yields a maximum 
value at some proper radius R. This maximum point is due to the strong boundary effect in the 
torsional resonator, as disclosed by the temperature distribution in the resonator. Resulting from 
the same boundary effect, the quality factor shows a power law temperature-dependence with power 
factors bellow 1.0. The theoretical results supply some valuable information for the manipulation 
of the quality factor in future experimental devices based on the torsional mechanical resonator. 
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I. INTRODUCTION 



Mechanical resonators have practical applications in 
fields hke accurate mass or force sensor The mechani- 
cal resonator is built based on a special natural vibration 
of the system. A gigahertz resonator was proposed the- 
oretically using multi-walled carbon nanotubes^ which 
was inspired by a retraction effect in this system.— This 
resonator is built on the oscillating vibration of multi- 
walled carbon nanotubes. The self-retracting effect, 
which is similar as resonator, was also observed in multi- 
ple graphene layers.— Using the cantilevered carbon nan- 
otube, several experimental groups have constructed a 
beam resonator, which actually takes advantage of the 
bending vibration (i.e. flexure mode).™!^!^, Besides the 
bending vibration, the torsional vibration of the single- 
walled carbon nanotube has also been used to create 
a torsional pendulum to measure the torsional force 
precisely.^ For graphene, both experimental^^— and the- 
oretical works! —1? have so far focused on the study of res- 
onators based on the bending vibration. However, there 
is still no study on possible mechanical resonator based 
on the torsional vibration of graphene. 



In this paper, we apply molecular dynamics (MD) sim- 
ulation to record the correlation function of the normal 
mode coordinate. The frequency uj and life time r of 
the torsional vibration are obtained from the correlation 
function. We then calculate the quality factor (Qf) of 
the torsional resonator hy Qp = ujt /2TT. Calculations 
focus on the size dependence and temperature depen- 
dence of the quality factor, which show the importance 
of the boundary effect. We also address an interesting re- 
current effect in finite sized resonators and estimate the 
corresponding recurrent speed. 
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FIG. 1: (Color online) Frequency versus the radius of the 
graphene torsional resonator. Insets are two resonator sam- 
ples with radius 4.5 A and 15.5 A. 



II. COMPUTATIONAL METHOD 

The structure of two resonators are shown as insets of 
Fig- HI The central six carbon atoms are fixed during the 
simulation. Free boundary condition is applied for car- 
bon atoms at the edge of the resonator. The arrow at- 
tached to each atom exhibits its vibrational amplitude in 
the torsional vibration. The maximum vibrational ampli- 
tudes are at the edge of the resonator. As a result, there is 
no stress at the edge, which is required by the free bound- 
ary condition. The torsional resonator can be actuated 
according to this vibrational morphology. The interac- 
tion between carbon atoms is described by the Brenner 
potential in the MD simulatioufii The Nose Hoover heat 
batbi^ii^ is applied to realize constant temperature when 
necessary. Newton equations are solved in terms of the 
velocity Verlet algorithm with a time step of 1 fs. 

A standard MD simulation for resonators is realized in 
three stagesJi--- In stage one, the resonator system is 
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thermalized to a constant temperature with the help of a 
heat bath. In stage two, the heat bath is removed and the 
mechanical resonator is actuated by generating the cor- 
responding natural vibration of the system. The phonon 
vibration is typically excited by displacing each atom ac- 
cording to the eigen vector of the vibration. In the final 
stage, the resonator system is allowed to oscillate freely, 
where the quality factor is calculated from the energy 
damping of the oscillation. We are not going to follow 
this standard procedure. Instead, we plan to perform MD 
simulation for the resonator in two steps. In step one, the 
system is thermalized at a constant temperature, which 
is realized by the Nose Hoover heat bath J^ii^ This step is 
the same as the first stage in the standard procedure. In 
step two, the heat bath is removed and the resonator is 
simulated in a microcanonical ensemble. In this step, we 
calculate the natural vibrations of the resonator system 
by solving an eigen value problem, and we can get the 
eigen vectors ^'^ for each phonon mode k. We then find 
the eigen vector corresponding to the torsional vibration 
^^torsiony Using this eigen vector, we can obtain 

the normal mode coordinate Q{t) of the torsional mode 
by following projection technique: 

3N 

where 3N is the total degree of freedom. {rj{t) — r^j) is the 
displacement of the j-th degree of freedom in the MD sim- 
ulation. Vj (t) is the position of the j-th degree of freedom 
at time t and is the equilibrium position. The time 
dependence of the normal mode coordinate, Q{t), is de- 
termined by the trajectory history from MD simulation. 
The quality factor of the torsional mechanical resonator 
is calculated by the definition Q p ~ vt — a;r/27rJi The 
frequency w and life time r of the torsional vibration can 
be derived from the correlation function. Under the sin- 
gle mode relaxation time approximation, the correlation 
function can be calculated by:^^ 

{Q%t)Q^Q))/{Q-{Q)Q%Q)) = cos{uj-t)e-'/-\ (1) 

where Q" (t) is the normal mode coordinate of the phonon 
mode a. Q'^iO) is the value at the beginning of the second 
step as described in above, uj'^ is the frequency and r'^ 
is the life time of this mode. There are two different 
interpretations for the average in the correlation function. 
It can be realized by the average over MD simulation 
time: 

M 

{Q^it„)Q^O)) = ^ E + tm)Q^tm), (2) 

m— 1 

where integer M is the total MD simulation step. M 
is usually in an order of 10'', so that the standard de- 
viation can be small enough. A long simulation time is 
required to reduce the statistical error due to the small 
size of the systems in our studied. Correlation functions 
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FIG. 2: (Color online) Correlation function of the torsional 
normal mode coordinate from averaging over simulation time. 
The blue curve is for a total simulation time of 10® simulation 
steps, while the black curve corresponds to a total simulation 
time of 10^ simulation steps. Panel (b) shows the correlation 
function in a larger time scale. 

for resonator with radius 4.5 A from this averaging with 
M — 10^ (in blue) and 10^ (in black) are shown in Fig. [51 
Panel (b) is for a larger time scale. This figure displays 
that the results show some dependence on the total sim- 
ulation time M, which may be due to possible statistical 
errors for large M. An interesting recurrent effect can be 
observed from panel (b). The correlation function does 
not decrease monotonously. Instead, it arises at some 
particular time tr after its first decaying. The recurrence 
is due to finite size effect and is not able to be eliminated 
by increasing simulation timeJ^ The physical picture of 
the recurrence is that the torsional phonon travels along 
the edge of the circular resonator and will come back to 
its starting point after some time. We will show further 
discussions on the recurrent effect in next section. 

Another interpretation for the average in the correla- 
tion function of Eq. ([l} is to average over initial condi- 
tions: 

1 ^ 

(O-(f„)g-(0)) = - ^ Qnt^)Qf{0), (3) 

i=l 

where i = 1,2, 3, N is the steady state from different 
initial conditions. This averaging method can help to 
avoid possible statistical error induced by extremely long 
simulation time. We will use this averaging approach in 
following. Correlation functions in resonator with radius 
4.5 A from this averaging method with N = 300 (in blue) 
and 450 (in black) are shown in Fig. |31 Panel (b) is for a 
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FIG. 3: (Color online) Correlation function of the torsional 
normal mode coordinate from averaging over initial condition 
of the molecular dynamics simulation. The blue curve is the 
average over 300 different initial conditions, while the black 
curve corresponds to an average over 450 different initial con- 
ditions. Panel (b) shows the correlation function in a larger 
time scale. 



larger time scale. The recurrent effect also exists in this 
averaging method as shown in panel (b). The frequency 
and life time are obtained by fitting Eq. ([ij to the corre- 
lation functions from MD simulation in Fig. [3] (a). The 
recurrent effect can be partially avoided by finishing the 
fitting procedure before the recurrence happens. That is 
why panel (a) is used in the fitting procedure rather than 
panel (b). In the following, the fitting procedure is done 
for the decay of the correlation function before recurrent 
time tr- Fig. m shows the frequency and life time for res- 
onator with radius 4.5 A from different sampling number 
A^. The value of these two quantities in the figure have 
been rescaled by the value at = 450. The frequency 
has almost the same value in the whole range of A. The 
life time reaches a saturated value after A^ > 300. We 
always use A^ = 350 for following calculations. 



III. RESULTS AND DISCUSSION 

Fig.[l]shows the natural frequency at K for resonators 
of different radii. The radius dependence of the frequency 
can be well described by function ut = 2628/i?^, where 
R is the radius of the resonator in A. This behavior is 
similar as a circular elastic solid thin plate under the 
same boundary condition. Treating the graphene res- 
onator as a circular elastic thin plate, we can obtain 

its torsional vibrational frequency: w = -mx/ e"^" '"^^ 
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FIG. 4: (Color online) The value of the frequency (uj) and life 
time (r) of the torsional resonator with radius 4.5 A versus 
the sampling number of the initial conditions. The displayed 
data is the original value rescaled by the value at sampling 
number 450. 



where m„ is the mass of neutron. 5*0 = 2.62A^ is 
the area of the two- atom unit cell in graphene sheet 
The effective torsional constant kg is estimated to be 
about 182.55 eV from the data in Fig. [TJ This value is 
larger than the Young's modules of the graphene which 
is about 50 eV^^"— indicating stronger stability for res- 
onator under torsional strain than the uniaxial strain. 
In usual, higher operational (resonant) frequency 
is favorable for the sensitivity of nanoscale res- 
onators. A sensitive mass sensor with resolution 
10 18 g can be obtained from carbon-nanotube- 
based electromechanical resonators with the op- 
erational frequency over 1.3 GHz.^'* Fig. [1] shows 
that the resonant frequency is over 1.0 GHz for 
torsional resonators with radius smaller than 70 
nm. This is helpful for the application of torsional 
resonators as sensitive mass sensor. 

Fig. O displays the life time at room temperature for 
the seven resonators with radius varying from 4.5 A to 
16 A. These points fall in a curve described by the func- 
tion T = l/(22i?"^ + 0.004) ps. The life time increases 
with increasing radius, which implies strong boundary 
scattering effect in the resonator. The boundary effect is 
of particular importance for free edges It is not easy to 
perform MD simulation for very large resonators, yet the 
fitting expression for the life time can be combined with 
the frequency w = 2628/i?^ to predict the quality factor 
of some large resonators at room temperature. This pre- 
diction is shown as inset (a) in the figure. The frequency 
formula at K has been used to predict the room temper- 
ature quality factor, since to is not sensitive to tempera- 
ture. The quality factor at 50 K has also been obtained 
by similar procedure. Both curves in inset (a) display 
a maximum value for the quality factor at some proper 
radius. For smaller radius, the boundary scattering is 
more important than the phonon-phonon scattering, so 



4 




radius (A) 



301.2 



^ 300.8 

&_ 

■(5 300.4 

&_ 

(1) 

Q. 

E 

(D 

300 



299.6 




4 8 12 

radial distance (A) 



16 



FIG. 5: (Color online) The room temperature life time versus 
the radius of the torsional resonator. Inset (a) is the quality 
factor from Qf ~ cjr/27r, with cj from the formula in Fig. [1] 
and life time by the fitting formula. Inset (b) is the difference 
between the recurrent time tr — 2-kR/0.5 and the life time r 
at room temperature. The x axis of both insets are the same 
as the main figure. 
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FIG. 6: (Color online) Temperature dependence for different 
physical quantities in the torsional resonator, (a). Recurrent 
velocity in resonator with radius 4.5 A. (b). The frequency 
in resonators with three different radius, (c). The life time of 
the torsional vibration mode in these three resonators, (d). 
The quality factor of these three torsional resonators. 



the increase of life time overpasses the decrease of fre- 
quency, leading to the increase of quality factor. In large 
resonator, the phonon-phonon scattering is more impor- 
tant, so the situation becomes opposite and the quality 
factor decreases with increasing radius. For a given ra- 
dius, the phonon-phonon scattering is stronger at higher 
temperature. Hence, the phonon-phonon scattering will 
become more important than the boundary scattering in 
resonators with smaller radii and at higher temperature. 
As a result, the maximum quality factor is achieved at 
smaller radius if temperature is higher. 

With the recurrent time tr and the circumference of 



FIG. 7: (Color online) The temperature distribution in the 
torsional resonator with radius 15.5 A. Inset shows the spatial 
temperature distribution by color. 



resonator L ~ we can estimate the recurrent ve- 

locity V = L/tr- Fig. [6] (a) shows the temperature de- 
pendence for the recurrent velocity of the resonator with 
R = 4.5 A. The recurrent velocity is larger at higher tem- 
perature, which is a common phenomenon in the nature. 
The curve can be fitted to a power function v = 0.09T°-^^ 
km/s. V ~ 0.5 km/s at room temperature. This value 
is close to the sound speed in the air, while it is much 
smaller than the longitudinal and transverse acoustic ve- 
locities in the graphene,— which is probably related to 
the edge effect in the torsional graphene resonator. 

We are interesting to investigate how the recurrence 
affects the correlation function, although we have only 
used the function for t < tr- In case of tr < r, the 
recurrent effect happens before the damping of the tor- 
sional vibration, so the recurrent wave will couple to the 
original wave in the correlation. This effect causes some 
unavoidable error for the correlation function after t > tr- 
That's another reason why we finish the fitting process 
before t < tr- The inset (b) in Fig. [5] shows the curve 
of tr — T at room temperature, where tr = 2i:R/v with 
constant v = 0.5 km/s as the recurrent velocity at 300 K. 
It shows that the recurrent effect will have influence for 
resonator with R < 180 A. If the MD simulation can be 
done for resonator with radius R > 180 A, the recurrent 
effect will be much weaker on the fitting procedure of 
the correlation function. For extremely large resonator, 
there will be no recurrent effect as tr is almost infinity. 

Fig. [5] (b) shows the temperature dependence for the 
frequency of three resonators with radius R = 4.5, 9.0, 
and 15.0 A. For all three resonators, a blue-shift of 
the frequency is observed with the increase of temper- 
ature. These three curves can be fitted to functions: 
uj/ujo = 1 + a* T^, with parameters (a, 6) as (0.01,0.5), 
(0.02,0.36), and (0.02,0.38). Overall, the frequency is 
not sensitive to temperature and the shift is within 20% 
below 300 K. We expect this blue-shift can be measured 
in the experiment, while current experiments focus on the 
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temperature dependence of other vibration modes such as 
radial breathing modes in the carbon nanotube^^^ Panel 
(c) is the life time at different temperatures. They can be 
described by power functions t — a*T~^, where parame- 
ters (a, b) are (23.8, 0.36), (186.3, 0.36), and (2255.8, 0.55) 
in three resonators. It is interesting that the power fac- 
tor b is less than 1.0 in all of the three resonators. The 
phonon life time is inversely proportional to the temper- 
ature in bulk materials?^ because the three-phonon scat- 
tering is the major process for the damping of the phonon 
vibration in bulk materials. However, in nano-materials 
such as the graphene resonator, the boundary effect pro- 
vides another channel for phonon decaying, which is tem- 
perature independent. The combination of the phonon- 
phonon scattering and the boundary effect leads to a 
power factor b less than 1.0. For larger resonator with 
R = 15.0 A, the power factor b = 0.55 is closer to the 
value of 1 than the other two smaller resonators where 
6 = 0.36. It implies that the boundary effect is weaker in 
larger resonator, which is quite reasonable. As a result of 
the temperature dependence of the life time, the quality 
factor shows similar behavior in panel (d) , where the two 
parameters (a, 6) in the power function are (92.7, 0.33), 
(178.9, 0.36), and (739.7, 0.53). Due to the boundary 
effect, the quality factor is very low for the torsional res- 
onator simulated in our work. The boundary effect on 
the quality factor has also been pointed out by Kim and 
Parkji^ If the atoms on the boundary are fixed and 
let the atoms at the center free, the resonator will 
undergo a bending oscillation instead of torsional 
oscillation. In this type of circular bending res- 
onator, Kim and Park have shown that the scat- 
tering by free edges will be considerably removed. 
As a result, the quality factor can be higher than 
2000 at room temperature and increase to be lO'^ 
at 10 In Fig. [71 the boundary effect is disclosed by 

the temperature distribution in the resonator. The tem- 
peratures of the inner six carbon atoms are zero, as they 
are fixed during MD simulation. There is a tempera- 
ture jump between the edge atoms and the inner carbon 
atoms, which results from the edge modes localized in 
this free boundary regioui^i^ 

We further remark that the torsional resonator de- 
scribed in our manuscript may be difficult to be realized 
in the experiment currently. The most difficult part is to 
fix the central region. A possible idea may be to fix the 
central region by hybridizing it to a silicon nanowire or 
diamond nanowire. Once the central region is fixed, it 
is straight forward to generate the torsional resonator by 
displacing the edge of the resonator using scanning elec- 
tron microscope, which has been successfully applied to 
push the surface graphene layer away from the graphitei^ 



Although the torsional resonator is not so easy to be real- 
ized with current experimental techniques, it is possible 
to be produced in the future with the development of the 
technique and some smart skills from experimentalists. 
Forgive the detailed technical difficulties in preparing the 
device, the significant contribution from our work is that 
we have proposed another theoretical approach to study 
the quality factor of the nanoresonator. Currently, there 
is only one MD simulation method for the investigation 
of mechanical resonators Our method can be used 
to obtain the pure intrinsic damping mechanism for the 
nanoresonator, as the quality factor is calculated from an 
equilibrium MD simulation and no mechanical oscillation 
is generated. 

It should be mentioned that we study the intrinsic 
damping mechanisms for the resonator, which is due to 
the phonon phonon scattering. This intrinsic damping 
mechanism depends on the structure and the interaction 
of the system. That is why we can calculate the quality 
factor from the correlation function of the normal mode 
coordinate. The phonon phonon scattering is described 
by the Brenner potential, which can successfully cap- 
ture both linear and nonlinear properties of the carbon- 
based system. 



IV. CONCLUSION 

To conclude, we perform MD simulation to investigate 
graphene-based torsional resonators. The quality factor 
is calculated hy Qp = a;r/27r, where w and r are the fre- 
quency and life time of the torsional vibration, uj and r 
are obtained from the correlation function of the normal 
mode coordinate. As a result of the boundary effect, 
the quality factor is small and has a maximum value 
for resonator with proper radius. The maximum qual- 
ity factor is achieved with smaller radius at higher tem- 
perature resulting from the interplay between boundary 
effect and the phonon-phonon scattering. The quality 
factor depends on temperature as a power law function 
Qp = a* with the power factor 6 < 1, which is due 
to the same boundary effect. We also investigate the re- 
current effect in finite sized mechanical resonators and 
the recurrent speed is found to be around 0.5 km/s at 
room temperature. 
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